home *** CD-ROM | disk | FTP | other *** search
/ ADA Programming Guide / ADA Programming Guide.iso / ada_gwu / arith.c < prev    next >
C/C++ Source or Header  |  1996-01-30  |  47KB  |  2,081 lines
  1. /*
  2.  * Copyright (C) 1985-1992  New York University
  3.  * 
  4.  * This file is part of the Ada/Ed-C system.  See the Ada/Ed README file for
  5.  * warranty (none) and distribution info and also the GNU General Public
  6.  * License for more details.
  7.  
  8.  */
  9.  
  10. /* Avoid problems with library routines on PC that require pointer
  11.  * to double instead of double (this is violation of ANSI specs).
  12.  * for now do not compile these modules
  13.  */
  14.  
  15. /* 
  16.  * 11-oct-85    shields
  17.  * fix so can handle most negative integer properly. The changes
  18.  * are more in the nature of a patch than a solution. The proper
  19.  * way is NOT to convert negative values to positive form before
  20.  * processing; instead, positive values should be converted to
  21.  * negative and conversions done on negative values. However, this
  22.  * can be put off until later.
  23.  *
  24.  * 5-jun-85    shields
  25.  * add rat_tol and int_tol which are like rat_toi and int_toi, resp.,
  26.  * except that they return long results. These are needed to support
  27.  * fixed-point in generator. 
  28.  *
  29.  * 1-aug-85    shields
  30.  * add calls to int_free() to free intermediate values known to be dead.
  31.  * add some copies where needed in building new rational values.
  32.  */
  33.  
  34. #include <stdlib.h>
  35. #include <stdio.h>
  36. #include <ctype.h>
  37. #include <string.h>
  38. #include <math.h>
  39. #include "config.h"
  40. #include "miscp.h"
  41.  
  42. typedef struct Rational_s {
  43.     int       *rnum;        /* numerator */
  44.     int       *rden;        /* denominator */
  45. }                Rational_s;
  46. typedef Rational_s * Rational;
  47.  
  48. #include "arithp.h"
  49.  
  50. static int *int_frl(long);
  51. static int *int_ptn(int);
  52. static int *subu_int(int *, int *);
  53. static void int_div(int *, int *, int **, int **);
  54. static int *int__addu(int *, int *);
  55. static int *int__norm(int *);
  56. static int *int_new(int);
  57. static void int_free(int *);
  58. static void rat_free(Rational);
  59. static double pow2(int);
  60.  
  61. /* Some of the procedures want to signal overflow by returning the
  62.  * string 'OVERFLOW'. In C we do this by setting global arith_overflow
  63.  * to non-zero if overflow occurs, zero otherwise
  64.  */
  65. int arith_overflow = 0;
  66.  
  67.  
  68. int    ADA_MIN_INTEGER;
  69. int    ADA_MAX_INTEGER;
  70. int    *ADA_MAX_INTEGER_MP;
  71. int    *ADA_MIN_INTEGER_MP;
  72. long    ADA_MIN_FIXED, ADA_MAX_FIXED;
  73. int    *ADA_MIN_FIXED_MP, *ADA_MAX_FIXED_MP;
  74. /* _LONG form is form that can be held in C (long) */
  75. #ifdef MAX_INTEGER_LONG
  76. long    MIN_INTEGER_LONG;
  77. int    *MAX_INTEGER_LONG_MP;
  78. int    *MIN_INTEGER_LONG_MP;
  79. #endif
  80.  
  81. #define ABS(x) ((x)<0 ? -(x) : (x))
  82. #define SIGN(x) ((x)<0 ? -1 : (x) == 0 ? 0 : 1)
  83.  
  84.  
  85. /*
  86.  * M U L T I P L E     P R E C I S I O N     I N T E G E R
  87.  *
  88.  *       A R I T H M E T I C       P A C K A G E
  89.  *
  90.  *             Robert B. K. Dewar
  91.  *              June 16th, 1980
  92.  *
  93.  * This package of routines implements multiple precision integer
  94.  * arithmetic using what are called the "classical algorithms" in
  95.  * Knuth "Art of Programming", Volume 2, Section 4.3.1. The style
  96.  * of the algorithms follows Knuth fairly closely, and section
  97.  * 4.3.1 can be consulted for further theoretical details.
  98.  *
  99.  * Multiple precision integers are represented as tuples of small
  100.  * integers in the range 0 to BAS - 1, where BAS is a power of 10
  101.  *(actually 10 ** DIGS) which is the base used to represent the
  102.  * long integers. Essentially the successive elements of the tuple
  103.  * are the digits of the representation in base BAS. All integers
  104.  * are normalized so that the first digit is non-zero, except in
  105.  * the case of zero itself. The sign is carried with the first
  106.  * digit value, all remaining digits are always positive.
  107.  *
  108.  * Some examples of the representation, assuming DIGS = 4 and
  109.  * BAS = 10000(note that the choice of BAS as a power of 10
  110.  * is for convenience of the conversion routines, and is not
  111.  * required for correct operation of the arithmetic algorithms).
  112.  
  113.  *     -123456789     [-1 , 2345 , 6789]
  114.  *        0     [0]
  115.  *    123456789     [1 , 2345, 6789]
  116.  
  117.  * The constants BAS and DIGS must be defined as global constants
  118.  * in a program using these routines. It is assumed that the value
  119.  * BAS * BAS - 1 can be represented as a SETL integer value.
  120.  
  121.  * The following routines are provided:
  122.  
  123.  *    int_abs        integer absolute value
  124.  *    int_add        integer addition
  125.  *    int_div        integer division
  126.  *    int_eql        integer equal to
  127.  *    int_exp        raise multiple integer to multiple integer power
  128.  *    int_fri        convert multiple precision integer from integer
  129.  *    int_frs        convert multiple precision integer from string
  130.  *    int_geq        integer greater than or equal to
  131.  *    int_gtr        integer greater than
  132.  *    int_len        number of digits in multiple precision integer
  133.  *    int_leq        integer less than or equal to
  134.  *    int_lss        integer less than
  135.  *    int_mod        integer modulus
  136.  *    int_mul        integer multiplication
  137.  *    int_neq        integer not equal to
  138.  *    int_ptn        integer power of ten
  139.  *    int_quo        integer quotient
  140.  *    int_rem        integer remainder
  141.  *    int_sub        integer subtraction
  142.  *    int_toi        convert integer to C integer
  143.  *    int_tos        integer to string
  144.  *    int_umin    integer unary minus
  145.  */
  146.  
  147. /* Internal procedures have names starting with int__ */
  148.  
  149. int    *int_abs(int *u)                                                /*;int_abs*/
  150. {
  151.     /* Absolute value of multiple precision integer */
  152.     int *pu;
  153.  
  154.     pu = int_copy(u);
  155.     if (pu[1] < 0)
  156.         pu[1] = -pu[1];
  157.     return pu;
  158. }
  159.  
  160. int    *int_add(int *u, int *v)                                        /*;int_add*/
  161. {
  162.     /* Add signed integers
  163.      * This procedure implements the familiar algorithm of comparing
  164.      * the signs, if the signs are the same, then the result is the
  165.      * sum of the magnitudes with the sign of the operands. If the
  166.      * signs differ, then the number with the smaller magnitude is
  167.      * subtracted from the larger magnitude and the result sign is
  168.      * that of the operand with the larger magnitude.
  169.      */
  170.  
  171.     int       *t;
  172.  
  173.     if (u[1] >= 0 && v[1] >= 0)
  174.         return (int__addu(u, v));
  175.     else
  176.         if (u[1] < 0 && v[1] < 0) {
  177.             u[1] = -u[1];
  178.             v[1] = -v[1];
  179.             t = int__addu(u, v);
  180.             u[1] = -u[1];
  181.             v[1] = -v[1];
  182.             t[1] = -t[1];
  183.             return t;
  184.         }
  185.         else {
  186.             int        us, vs;
  187.             if (us = (u[1] < 0)) {
  188.                 u[1] = -u[1];
  189.             }
  190.             if (vs = (v[1] < 0)) {
  191.                 v[1] = -v[1];
  192.             }
  193.             if (int_gtr(u, v)) {
  194.                 t = subu_int(u, v);
  195.                 if (vs) v[1] = -v[1];
  196.                 if (us) {
  197.                     u[1] = -u[1];
  198.                     t[1] = -t[1];
  199.                 }
  200.                 return t;
  201.             }
  202.             else {
  203.                 t = subu_int(v, u);
  204.                 if (us) u[1] = -u[1];
  205.                 if (vs) {
  206.                     v[1] = -v[1];
  207.                     t[1] = -t[1];
  208.                 }
  209.                 return t;
  210.             }
  211.         }
  212. }
  213.  
  214. int    int_eql(int *u, int *v)                                        /*;int_eql*/
  215. {
  216.     /* Compare multiple precision integers for equality */
  217.  
  218.     int        n;
  219.  
  220.     if (u[0] != v[0])
  221.         return FALSE;
  222.  
  223.     for (n = u[0]; n > 0; n--)
  224.         if (u[n] != v[n]) return FALSE;
  225.     return TRUE;
  226. }
  227.  
  228. int    *int_exp(int *u, int *v)                                        /*;int_exp*/
  229. {
  230.     /* Raise a multiple precision integer to a multiple precision integer
  231.      * power using a modified version of the Russian Peasant algorithm
  232.      * for exponentiation.
  233.      */
  234.  
  235.     int        sn;
  236.     int        u1;
  237.     int        i;
  238.     int        carry;
  239.     int        half;
  240.     int       *running, *runningp;
  241.     int       *result, *resultp;
  242.  
  243.     /* Compute sign of result: positive if v is even, the sign of u if
  244.      * v is odd.
  245.      */
  246.  
  247.     sn =((v[v[0]] % 2) == 1) ? SIGN(u[1]) : 1;
  248.     u1 = u[1];
  249.     if (u[1] < 0)
  250.         u[1] = -u1;
  251.     v = int_copy(v);
  252.  
  253.     /* Starting the result at 1 and running at u, loop through the binary
  254.      * digits of v, squaring running each time, and multiplying the result
  255.      * by the current value of running each time a 1-bit is found.
  256.      */
  257.  
  258.     result = int_con(1);
  259.     running = int_copy(u);
  260.  
  261.     while(int_nez(v)) {
  262.         /* If v is odd then multiply result by running. */
  263.  
  264.         if (v[v[0]] % 2 == 1) {
  265.             resultp = result; /* save current value */
  266.             result = int_mul(result, running);
  267.             int_free(resultp); /* free dead value */
  268.         }
  269.  
  270.         /* Square running. */
  271.  
  272.         runningp = running; /* save current value */
  273.         running = int_mul(running, running);
  274.         int_free(runningp); /* free dead value */
  275.  
  276.         /* Halve v. */
  277.  
  278.         carry = 0;
  279.         for (i = 1; i <= v[0]; i++) {
  280.             half = BAS * carry + v[i];
  281.             carry = half % 2;
  282.             v[i] = half / 2;
  283.         }
  284.         v = int__norm(v);
  285.     }
  286.  
  287.     int_free(running);
  288.     int_free(v);
  289.  
  290.     if (sn < 1)
  291.         result[1] = -result[1];
  292.     u[1] = u1;
  293.     return result;
  294. }
  295.  
  296. int    *int_fri(int i)                                                /*;int_fri*/
  297. {
  298.     /* Convert an integer to a multiple precision integer.
  299.      *
  300.      * First check the sign of i.
  301.      */
  302.  
  303.     int        sn = 0;
  304.     int        n = i;
  305.     int       *t;
  306.     int        ti;
  307.     int        d;
  308.  
  309.     if (i < 0) {
  310.         /* Special test for most negative as code below won't work
  311.          * for this single value (on twos-complement machines)
  312.          */
  313.         if (i == ADA_MIN_INTEGER)
  314.             return int_copy(ADA_MIN_INTEGER_MP);
  315.         sn = -1;
  316.         n = -i;
  317.     }
  318.     /* Determine length of result */
  319.     d = 0;
  320.     while(n) {
  321.         d++;
  322.         n /= BAS;
  323.     }
  324.     /* Now build up t in groups of BAS digits until i is depleted. */
  325.  
  326.     if (i < 0)
  327.         i = -i;
  328.     if (i > 0) {
  329.         t = int_new(d);
  330.         ti = t[0];
  331.         while(i) {
  332.             t[ti--] = i % BAS;
  333.             i /= BAS;
  334.         }
  335.     }
  336.     else
  337.         t = int_con(0);
  338.  
  339.     if (sn < 0)
  340.         t[1] = -t[1];
  341.  
  342.     return t;
  343. }
  344.  
  345. #ifdef MAX_INTEGER_LONG
  346. /* variant of int_fri needed for PC only */
  347. static int *int_frl(long i)                                        /*;int_frl*/
  348. {
  349.     /* Convert a long integer to a multiple precision integer.
  350.      *
  351.      * First check the sign of i.
  352.      */
  353.  
  354.     int        sn = 0;
  355.     long    n = i;
  356.     int       *t;
  357.     int        ti;
  358.     int        d;
  359.  
  360.     if (i < 0) {
  361.         /* Special test for most negative as code below won't work
  362.          * for this single value (on twos-complement machines)
  363.          */
  364.         sn = -1;
  365.         n = -i;
  366.     }
  367.     /* Determine length of result */
  368.     d = 0;
  369.     while(n) {
  370.         d++;
  371.         n /= BAS;
  372.     }
  373.     /* Now build up t in groups of BAS digits until i is depleted. */
  374.  
  375.     if (i < 0)
  376.         i = -i;
  377.     if (i > 0) {
  378.         t = int_new(d);
  379.         ti = t[0];
  380.         while(i) {
  381.             t[ti--] = i % BAS;
  382.             i /= BAS;
  383.         }
  384.     }
  385.     else
  386.         t = int_con(0);
  387.  
  388.     if (sn < 0)
  389.         t[1] = -t[1];
  390.  
  391.     return t;
  392. }
  393. #else
  394. /* if ints and longs same size, just use int_fri */
  395. static int    *int_frl(long i)                                        /*;int_frl*/
  396. {
  397.     return int_fri((int)i);
  398. }
  399. #endif
  400.  
  401. int    *int_frs(char *s)                                            /*;int_frs*/
  402. {
  403.     /* Convert a string to multiple precision integer format.  The string
  404.      * s is a non-empty sequence of digits, possibly preceded by a sign
  405.      *(+ or -).
  406.      *
  407.      * Erroneous strings are converted to OM.
  408.      *
  409.      * Since the base is a power of ten, the process of conversion
  410.      * amounts simply to converting groups of DIGS digits.
  411.      */
  412.  
  413.     char   *t;
  414.     int        dg;
  415.     int       *r, len;
  416.     int        ri, z, sn;
  417.     int        n;
  418.  
  419.     if (*s == '+') {
  420.         sn = 0;
  421.         s++;
  422.     }
  423.     else
  424.         if (*s == '-') {
  425.             sn = -1;
  426.             s++;
  427.         }
  428.         else
  429.             sn = 0;
  430.  
  431.     /* now determine length */
  432.     t = s;
  433.     dg = 0;
  434.     while(*t) {
  435.         if (!isdigit(*t++))
  436.             return (int *)0;
  437.         dg++;
  438.     }
  439.     if (dg == 0)
  440.         return (int *)0;
  441.  
  442.     z = dg % DIGS;
  443.     r = int_new((dg / DIGS) +(z ? 1 : 0));
  444.     ri = 1;
  445.     len = z ? z : DIGS;
  446.  
  447.     while(dg > 0) {
  448.         n = 0;
  449.         dg -= len;
  450.         while(len--) {        /* convert next len digits */
  451.             n = n * 10 +(*s++ - '0');
  452.         }
  453.         len = DIGS;
  454.         r[ri++] = n;
  455.     }
  456.     if (sn < 0)
  457.         r[1] = -r[1];
  458.     return int__norm(r);
  459. }
  460.  
  461. int    int_geq(int *u, int *v)                                        /*;int_geq*/
  462. {
  463.     /* Compare multiple precision integers: return true if u >= v, false
  464.      * otherwise.
  465.      */
  466.  
  467.     return int_eql(u, v) || int_gtr(u, v);
  468. }
  469.  
  470. int    int_gtr(int *u, int *v)                                        /*;int_gtr*/
  471. {
  472.     /* Compare multiple precision integers: return true if u > v, false
  473.      * otherwise.
  474.      */
  475.  
  476.     int        i, neg;
  477.  
  478.     /* signs are different                       */
  479.  
  480.     if (u[1] >= 0 && v[1] < 0) return TRUE;
  481.  
  482.     if (u[1] < 0 && v[1] >= 0) return FALSE;
  483.  
  484.     /* Now we have a real compare(signs the same) */
  485.  
  486.     neg = u[1] < 0;   /* get the sign */
  487.  
  488.     if (u[0] > v[0]) return (neg ? FALSE : TRUE);
  489.  
  490.     if (u[0] < v[0]) return (neg ? TRUE : FALSE);
  491.  
  492.     /* Now the lengths are the same  */
  493.  
  494.     if(u[1] != v[1]) return (u[1] > v[1]);
  495.  
  496.     /* Most significant digit is the same */
  497.  
  498.     for (i = 2; i <= u[0]; i++) {
  499.         if (u[i] != v[i]) {
  500.             return (neg ?(v[i] > u[i]) :(u[i] > v[i]));
  501.         }
  502.     }
  503.  
  504.     /* Numbers are the same */
  505.  
  506.     return FALSE;
  507. }
  508.  
  509. /* Not used */
  510. int    int_len(int *u)                                                /*;int_len*/
  511. {
  512.     /* Return the number of digits in a multiple precision integer. */
  513.  
  514.     int        n = 1;
  515.     int        v;
  516.  
  517.     v = u[1];
  518.     if (v < 0)
  519.         v = -v;            /* take absolute value */
  520.     while(v) {
  521.         n++;
  522.         v /= 10;
  523.     }
  524.     return n +(u[0] - 1) * DIGS;
  525. }
  526.  
  527. /* Not used */
  528. int int_leq(int *u, int *v)                                            /*;int_leq*/
  529. {
  530.     /* Compare multiple precision integers: return true if u <= v, false
  531.      * otherwise.
  532.      */
  533.  
  534.     return ! int_gtr(u, v);
  535. }
  536.  
  537. int    int_lss(int *u, int *v)                                        /*;int_lss*/
  538. {
  539.     /* Compare multiple precision integers: return true if u < v, false
  540.      * otherwise.
  541.      */
  542.  
  543.     return ! int_geq(u, v);
  544. }
  545.  
  546. int    *int_mod(int *u, int *v)                                        /*;int_mod*/
  547. {
  548.     /* Find u mod v where the mod operation is defined as:
  549.      *
  550.      *    u mod v = u - v * N    such that sign(u mod v) = sign v
  551.      *                      and  abs(u mod v) < abs v
  552.      */
  553.  
  554.     int       *r, *s;
  555.  
  556.     r = int_rem(u, v);
  557.  
  558.     if (r == (int *)0 || r[1] == 0 ||(SIGN(u[1]) == SIGN(v[1])))
  559.         return r;
  560.     else {
  561.         s = int_add(r, v);
  562.         int_free(r);
  563.         return s;
  564.     }
  565. }
  566.  
  567. int    *int_mul(int *u, int *v)        /*;int_mul*/
  568. {
  569.     /* Multiply signed integers, cf Knuth 4.3.1 Algorithm M
  570.      *
  571.      * Multiplication is similar to, but not identical with, the
  572.      * usual pencil and paper algorithm. The main difference is
  573.      * that the sum is accumulated as we go along, rather than
  574.      * forming all the partial sums and adding them up at the end.
  575.      *
  576.      * First we acquire the sign of the result, and set each number
  577.      * to its absolute value, thus the multiplication proper always
  578.      * works with non-negative integers.
  579.      */
  580.  
  581.     int        sn;
  582.     int        u1;
  583.     int        v1;
  584.     int       *w;
  585.     int        i, j;
  586.     int        k;
  587.     int        t;
  588.  
  589.     sn = u[1] * v[1];
  590.     u1 = u[1];
  591.     v1 = v[1];
  592.     if (u[1] < 0)
  593.         u[1] = -u1;
  594.     if (v[1] < 0)
  595.         v[1] = -v1;
  596.     /* Initialize result to all zeroes(actually it is only absolutely
  597.      * necessary to initialize the last #v digits of w to zero).
  598.      */
  599.  
  600.     w = int_new(u[0] + v[0]);
  601.  
  602.     /* Outer loop, through digits of multiplier */
  603.  
  604.     for (j = v[0]; j > 0; j--) {
  605.         /* The inner loop, through the digits of the multiplicand, is
  606.          * similar to the inner loop of the addition, except that the
  607.          * product is added in, and the carry, k, can exceed 1.
  608.          */
  609.  
  610.         k = 0;
  611.         for (i = u[0]; i > 0; i--) {
  612.             t = u[i] * v[j] + w[i + j] + k;
  613.             w[i + j] = t % BAS;
  614.             k = t / BAS;
  615.         }
  616.  
  617.         /* The final step in the inner loop is to store the final
  618.          * carry in the next position in the result.
  619.          */
  620.  
  621.         w[j] = k;
  622.     }
  623.  
  624.     /* The result must be normalized(there could be one leading zero),
  625.      * and then the result sign attached to the returned value.
  626.      */
  627.  
  628.     w = int__norm(w);
  629.     /* Restore arguments to entry value */
  630.     u[1] = u1;
  631.     v[1] = v1;
  632.     if (sn < 0)
  633.         w[1] = -w[1];
  634.     return w;
  635. }
  636.  
  637. /* Not used */
  638. int    int_neq(int *u, int *v)                                        /*;int_neq*/
  639. {
  640.     /* Compare multiple precision integers for inequality */
  641.  
  642.     return ! int_eql(u, v);
  643. }
  644.  
  645. static int    *int_ptn(int n)                                        /*;int_ptn*/
  646. {
  647.     /* Return the result 10**n as a multiple precision
  648.      * integer, where n is a non-negative simple integer.
  649.      */
  650.  
  651.     int        d1;
  652.     int       *p;
  653.     int        i;
  654.  
  655.     d1 = 1;
  656.     for (i = 1; i <= (n % DIGS); i++)
  657.         d1 *= 10;
  658.     p = int_new(1 +(n / DIGS));
  659.     p[1] = d1;
  660.     return p;
  661. }
  662.  
  663. int    *int_quo(int *u, int *v)                                        /*;int_quo*/
  664. {
  665.     /* Obtain quotient of dividing u by v */
  666.  
  667.     int       *q, *r;
  668.  
  669.     if (int_eqz(v))
  670.         return (int *)0;
  671.     int_div(u, v, &q, &r);
  672.     if(r != (int *)0) int_free(r); /* remainder not needed, free it */
  673.     return q;
  674. }
  675.  
  676. int    *int_rem(int *u, int *v)                                        /*;int_rem*/
  677. {
  678.     /* Obtain remainder from dividing u by v, where u rem v is defined as:
  679.      *
  680.      *    u rem v = u -(u/v) * v       such that sign(u rem v) = sign u
  681.      *                      and  abs(u rem v) < abs v
  682.      */
  683.  
  684.     int       *q, *r;
  685.  
  686.     if (int_eqz(v))
  687.         return (int *)0;
  688.     int_div(u, v, &q, &r);
  689.     if (q != (int *) 0) int_free(q); /* quotient not needed, free it */
  690.     return r;
  691. }
  692.  
  693. int    *int_sub(int *u, int *v)                                        /*;int_sub*/
  694. {
  695.     /* Subtract signed integers
  696.      *
  697.      * There is no point in duplicating the int_add code, so we
  698.      * simply negate the right argument and call the add routine!
  699.      */
  700.  
  701.     int *w;
  702.  
  703.     v[1] = -v[1];
  704.     w = int_add(u, v);
  705.     v[1] = -v[1];
  706.     return w;
  707. }
  708.  
  709. int    int_toi(int *t)                                                /*;int_toi*/
  710. {
  711.     /* Convert a multiple precision integer to a SETL integer, if possible.
  712.      * Otherwise, return 'OVERFLOW'.
  713.      *
  714.      * First check if it overflows.
  715.      */
  716.  
  717.     int        sgn;
  718.  
  719.     int        x;
  720.     int        i;
  721.  
  722.     arith_overflow = 0;        /* reset overflow flag */
  723.     sgn = SIGN(t[1]);
  724.  
  725.     if (sgn == 0)
  726.         return 0;
  727.     if (sgn > 0)
  728.         t[1] = -t[1];
  729.  
  730.     /* the value of t must be in the interval */
  731.     if (int_lss(t, ADA_MIN_INTEGER_MP)
  732.       || int_gtr(t, ADA_MAX_INTEGER_MP)) {
  733.         if (sgn > 0)
  734.             t[1] = -t[1];    /* restore first component */
  735.         arith_overflow = 1;
  736.         return 0;
  737.     }
  738.  
  739.     x = t[1]; /* set first part of result (must do here since negative) */
  740.     for (i = 2; i <= t[0]; i++)
  741.         x = BAS * x - t[i];
  742.  
  743.     if (sgn > 0) {
  744.         t[1] = -t[1];        /* restore sign if negative */
  745.         x = -x;            /* and give result proper value */
  746.     }
  747.  
  748.     return x;
  749. }
  750.  
  751. #ifdef MAX_INTEGER_LONG
  752. long int_tol(int *t)                                            /*;int_tol*/
  753. {
  754.     /* Convert a multiple precision integer to a long integer, if possible.
  755.      * Otherwise, return 'OVERFLOW'.
  756.      *
  757.      * First check if it overflows.
  758.      */
  759.  
  760.     int        sgn;
  761.     long        x;
  762.     long    res;
  763.     int        i;
  764.  
  765.     arith_overflow = 0;        /* reset overflow flag */
  766.     sgn = SIGN(t[1]);
  767.  
  768.     if (sgn == 0)
  769.         return (long)0;
  770.     if (sgn < 0) {
  771. #ifdef MAX_INTEGER_LONG
  772.         if (int_eql(t, MIN_INTEGER_LONG_MP)) 
  773.             return MIN_INTEGER_LONG;
  774. #else
  775.         if (int_eql(t, ADA_MIN_INTEGER_MP))
  776.             return ADA_MIN_INTEGER;
  777. #endif
  778.  
  779.         /* since not most negative, can change sign */
  780.         t[1] = -t[1];
  781.     }
  782.  
  783. #ifdef MAX_INTEGER_LONG
  784.     if (int_gtr(t, MAX_INTEGER_LONG_MP)) {
  785.         arith_overflow = 1;
  786.         return (long) 0;
  787.     }
  788. #else
  789.     if (int_gtr(t, ADA_MAX_INTEGER_MP)) {
  790.         arith_overflow = 1;
  791.         return (long) 0;
  792.     }
  793. #endif
  794.  
  795.     x = 0;
  796.     for (i = 1; i <= t[0]; i++)
  797.         x = BAS * x + t[i];
  798.  
  799.     if (sgn < 0)
  800.         t[1] = -t[1];        /* restore sign if negative */
  801.  
  802.     return (sgn * x);
  803. }
  804. #else
  805. /* if longs and ints are same size, no need for separate procedure */
  806. long int_tol(int *t)                                            /*;int_tol*/
  807. {
  808.     return (long) int_toi(t);
  809. }
  810. #endif
  811.  
  812. char   *int_tos(int *u)                                                /*;int_tos*/
  813. {
  814.     /* Convert a multiple precision integer to a string.
  815.      * The string is a non-empty digit sequence with a possible leading
  816.      * minus sign(but a plus sign is never generated).
  817.      *
  818.      * As in int_frs, the fact that the base is a power of ten means
  819.      * that the conversion is simply a matter of converting successive
  820.      * digits to decimal strings of length DIGS.
  821.      */
  822.  
  823.     char   *s, *t;
  824.     int        i, n;
  825.     if (u == (int *)0) chaos("int_tos: arg null *");
  826.  
  827.     n = u[0] * DIGS;
  828.     if (u[1] < 0)
  829.         n += 1;            /* if need minus sign */
  830.     s = t = emalloct((unsigned) n + 1, "int-to-s");
  831.     sprintf(s, "%d", u[1]);
  832.     for (i = 2; i <= u[0]; i++) {
  833.         while(*++s);        /* move to end of converted string */
  834.         sprintf(s, "%0*d", DIGS, u[i]);
  835.     }
  836.  
  837.     return t;
  838. }
  839.  
  840. int    *int_umin(int *u)                                        /*;int_umin*/
  841. {
  842.     /* Unary minus on multiple precision integer. */
  843.  
  844.     u = int_copy(u);
  845.     u[1] = -u[1];
  846.     return u;
  847. }
  848.  
  849. static int    *subu_int(int *u, int *v)                                /*subu_int*/
  850. {
  851.     /* Subtract unsigned integers, u >= v, cf Knuth 4.3.1 Algorithm S
  852.      *
  853.      *(note that we know as an entry condition that #u >= #v).
  854.      */
  855.  
  856.     int        ui, vi;
  857.     int       *w;
  858.     int        k;            /* borrow */
  859.  
  860.     w = int_new(u[0]);
  861.     ui = u[0];
  862.     vi = v[0];
  863.  
  864.     /* The subtraction is similar to the addition, except that k now
  865.      * represents the borrow condition and has values 0 or -1.
  866.      */
  867.  
  868.     k = 0;
  869.     while(vi) {
  870.         w[ui] = u[ui] - v[vi--] + k;
  871.         if (w[ui] < 0) {
  872.             w[ui] += BAS;
  873.             k = -1;
  874.         }
  875.         else
  876.             k = 0;
  877.         ui--;
  878.     }
  879.  
  880.     /* Loop over remaining digits(if any) of u */
  881.  
  882.     while(ui) {
  883.         w[ui] = u[ui] + k;
  884.         if (w[ui] < 0) {
  885.             w[ui] += BAS;
  886.             k = -1;
  887.         }
  888.         else
  889.             k = 0;
  890.         ui--;
  891.     }
  892.  
  893.     /* We cannot have a final borrow(since the entry condition
  894.      * required that u >= v). However, we must normalize the
  895.      * result since it is possible for leading zeroes to appear.
  896.      */
  897.  
  898.     w = int__norm(w);
  899.     return w;
  900. }
  901.  
  902. static void int_div(int *u, int *v, int **qa, int **ra)                /*;int_div*/
  903. {
  904.     /* Obtain quotient and remainder of signed integers,
  905.      * cf Knuth 4.3.1 Algorithm D
  906.      * Result is returned as a tuple [quotient, remainder].
  907.      *
  908.      * This is by far the most difficult of the four basic operations.
  909.      * This is because the paper and pencil algorithm involves certain
  910.      * amounts of guess work which cannot be programmed directly. The
  911.      * approach(which is analyzed in detail in Knuth) is to reduce
  912.      * the guess work by computing a rather good guess at each digit
  913.      * of the result, and then correcting if the guess is wrong.
  914.      *
  915.      * If the divisor is zero, return om.
  916.      */
  917.  
  918.     int        i, j, k, du, v1, p, c, d;
  919.     int        rr, rs;
  920.     int       *r, *q;
  921.     int       *t, *tt;
  922.     int        ti, n, m, qe, qs;
  923.  
  924.     if (int_eqz(v)) {
  925.         *qa = (int *)0;
  926.         *ra = (int *)0;
  927.         return;
  928.     }
  929.  
  930.     /* A special case, if u is shorter than v, the result is 0 */
  931.  
  932.     if (u[0] < v[0]) {
  933.         *qa = int_con(0);
  934.         *ra = int_copy(u);
  935.         return;
  936.     }
  937.  
  938.     /* Otherwise we initialize as in multiplication by obtaining the
  939.      * result sign and then working with non-negative integers.
  940.      */
  941.  
  942.     u = int_copy(u); /* arguments used destructively, so copy */
  943.     v = int_copy(v);
  944.     qs = u[1] * v[1];
  945.     rs = u[1];
  946.     v1 = v[1];
  947.     if (rs < 0)
  948.         u[1] = -rs;
  949.     if (v1 < 0)
  950.         v[1] = -v1;
  951.  
  952.     /* The case of a one digit divisor is treated specially. Not only
  953.      * is this more efficient, but the general algorithm assumes that
  954.      * the divisor contains at least two digits.
  955.      */
  956.  
  957.     if (v[0] == 1) {
  958.         q = int_new(u[0]);
  959.  
  960.         /* Basically this case is straight-forward. Since we can represent
  961.          * numbers up to BAS * BAS - 1, we can do the steps of the division
  962.          * exactly without any need for guess work. The division is then
  963.          * done left to right(essentially it is the short division case).
  964.          */
  965.         rr = 0;
  966.         for (j = 1; j <= u[0]; j++) {
  967.             du = rr * BAS + u[j];
  968.             q[j] = du / v[1];
  969.             rr = du % v[1];
  970.         }
  971.  
  972.         r = int_con(rr);
  973.     }
  974.     /* Here is where we must do the full long division algorithm */
  975.  
  976.     else {
  977.         n = v[0];
  978.         m = u[0] - v[0];
  979.         q = int_new(m+1);
  980.         t = int_new(u[0] + 1);    /* extend u */
  981.         for (j = 1; j <= u[0]; j++)
  982.             t[j + 1] = u[j];
  983.         int_free(u);
  984.         u = t;
  985.  
  986.         /* The first step is to multiply both the divisor and dividend
  987.          * by a scale factor. Obviously such scaling does not affect
  988.          * the quotient. The purpose of this scaling is to ensure that
  989.          * the first digit of the divisor is at least BAS / 2. This
  990.          * condition is required for proper operation of the quotient
  991.          * estimation algorithm used in the division loop. Note that
  992.          * we added an extra digit at the front of the dividend above.
  993.          */
  994.  
  995.         d = BAS /(v[1] + 1);
  996.  
  997.         c = 0;
  998.         for (j = u[0]; j > 0; j--) {
  999.             p = u[j] * d + c;
  1000.             u[j] = p % BAS;
  1001.             c = p / BAS;
  1002.         }
  1003.  
  1004.         c = 0;
  1005.         for (j = v[0]; j > 0; j--) {
  1006.             p = v[j] * d + c;
  1007.             v[j] = p % BAS;
  1008.             c = p / BAS;
  1009.         }
  1010.  
  1011.         /* This is the major loop, corresponding to long division steps */
  1012.  
  1013.         for (j = 1; j <= (m + 1); j++) {
  1014.             /* Guess the next quotient digit by doing a division based on the
  1015.              * leading digits. This estimate is never low and at most 2 high.
  1016.              */
  1017.             qe =(u[j] != v[1])
  1018.                 ?((u[j] * BAS + u[j + 1]) / v[1])
  1019.                 :(BAS - 1);
  1020.  
  1021.             /* The following loop refines this guess so that it is almost always
  1022.              * correct and is at worst one too high(see Knuth for proofs).
  1023.              */
  1024.  
  1025.             while((v[2] * qe) >
  1026.                 ((u[j] * BAS + u[j + 1] - qe * v[1]) * BAS + u[j + 2])) {
  1027.                 qe -= 1;
  1028.             }
  1029.  
  1030.             /* Now(for the moment accepting the estimate as correct), we
  1031.              * subtract the appropriate multiple of the divisor. This is
  1032.              * similar to the inner loop of the multiplication routine.
  1033.              */
  1034.             c = 0;
  1035.             for (k = n; k > 0; k--) {
  1036.                 du = u[j + k] - qe * v[k] + c;
  1037.                 u[j + k] = du % BAS;
  1038.                 c = du / BAS;
  1039.                 if (u[j + k] < 0) {
  1040.                     u[j + k] += BAS;
  1041.                     c -= 1;
  1042.                 }
  1043.             }
  1044.             u[j] += c;
  1045.  
  1046.             /* If the estimate was one off, then u(j) went negative when the
  1047.              * final carry was added above. In this case, we add back the
  1048.              * divisor once, and adjust the quotient digit.
  1049.              */
  1050.  
  1051.             if (u[j] < 0) {
  1052.                 qe -= 1;
  1053.  
  1054.                 c = 0;
  1055.                 for (k = n; k > 0; k--) {
  1056.                     u[j + k] += v[k] + c;
  1057.                     if (u[j + k] >= BAS) {
  1058.                         u[j + k] -= BAS;
  1059.                         c = 1;
  1060.                     }
  1061.                     else
  1062.                         c = 0;
  1063.                 }
  1064.  
  1065.                 u[j] += c;
  1066.             }
  1067.  
  1068.             /* Store the next quotient digit */
  1069.  
  1070.             q[j] = qe;
  1071.         }
  1072.  
  1073.         /* Remainder is in u(m+2..m+n+1) but must be divided by d. */
  1074.  
  1075.  
  1076.         /* SETL:      r:=int_quo(u(m+2..m+n+1), [d]); */
  1077.         t = int_new(n);
  1078.         ti = 1;
  1079.         for (i = m + 2; i <= (m + n + 1); i++)
  1080.             t[ti++] = u[i];
  1081.         r = int_quo(t, tt=int_con(d));
  1082.         int_free(t);
  1083.         int_free(tt);
  1084.     }
  1085.  
  1086.     /* All done, except for normalizing the quotient
  1087.      * to remove leading zeroes and supplying the
  1088.      * proper result sign to the returned values.
  1089.      */
  1090.  
  1091.     q = int__norm(q);
  1092.     if (qs < 0)
  1093.         q[1] = -q[1];
  1094.     if (rs < 0)
  1095.         r[1] = -r[1];
  1096.     *qa = q;
  1097.     *ra = r;
  1098.  
  1099.     int_free(u);
  1100.     int_free(v);
  1101. }
  1102.  
  1103. static int    *int__addu(int *u, int *v)                        /*;int__addu*/
  1104. {
  1105.     /* Add unsigned integers, cf Knuth 4.3.1 Algorithm A
  1106.      *
  1107.      * We first do the addition just to see if final carry, so we
  1108.      * can determine correct length of result. Then we allocate
  1109.      * the result and perform the addition.
  1110.      */
  1111.  
  1112.     int        k = 0;        /* carry */
  1113.     int       *w;
  1114.     int        r;
  1115.     int        ui, vi, wi;
  1116.  
  1117.     /* Adjust so u points to longer argument */
  1118.  
  1119.     if (v[0] > u[0]) {        /* if second longer */
  1120.         w = u;
  1121.         u = v;
  1122.         v = w;
  1123.     }
  1124.     ui = u[0];            /* length of first argument */
  1125.     vi = v[0];
  1126.     while(vi) {
  1127.         r = u[ui--] + v[vi--] + k;
  1128.         if (r >= BAS)
  1129.             k = 1;
  1130.         if (r < BAS)
  1131.             k = 0;
  1132.     }
  1133.     while(ui) {
  1134.         r = u[ui--] + k;
  1135.         if (r >= BAS)
  1136.             k = 1;
  1137.         if (r < BAS)
  1138.             k = 0;
  1139.     }
  1140.     /* if k nonzero, result needs extra component */
  1141.     w = int_new(u[0] + k);
  1142.     /* Now perform addition for real */
  1143.  
  1144.     /* Now we perform the addition from right to left in the familiar
  1145.      * paper and pencil manner. The variable k is the carry from the
  1146.      * previous column, which has either the value zero or one.
  1147.      */
  1148.  
  1149.     ui = u[0];
  1150.     vi = v[0];
  1151.     wi = w[0];
  1152.     k = 0;
  1153.     while(vi) {
  1154.         w[wi] = u[ui--] + v[vi--] + k;
  1155.         if (w[wi] >= BAS) {
  1156.             w[wi] -= BAS;
  1157.             k = 1;
  1158.         }
  1159.         else
  1160.             k = 0;
  1161.         wi--;
  1162.     }
  1163.     /* Now add in remainder of longer number */
  1164.     while(ui) {
  1165.         w[wi] += u[ui--] + k;
  1166.         if (w[wi] >= BAS) {
  1167.             w[wi] -= BAS;
  1168.             k = 1;
  1169.         }
  1170.         else
  1171.             k = 0;
  1172.         wi--;
  1173.     }
  1174.     if (k)
  1175.         w[1] = 1;        /* if final carry */
  1176.  
  1177.     return w;
  1178. }
  1179.  
  1180. /* Note that int__norm does not alter its argument, but returns
  1181.  * pointer to(possibly normalized) multi-precision integer.
  1182.  */
  1183. static int    *int__norm(int *u)                                /*;int__norm*/
  1184. {
  1185.     /* Remove leading zeroes from calculated value
  1186.      *
  1187.      * The representation used in this package requires that all integer
  1188.      * values be normalized, i.e. the first digit of any stored value
  1189.      * must be non-zero except for the special case of zero itself. The
  1190.      * normalize routine is called to ensure this condition is met.
  1191.      *
  1192.      */
  1193.  
  1194.     int        i, j, *v;
  1195.  
  1196.     if (u[0] == 1 || u[1] != 0)
  1197.         return (u);
  1198.  
  1199.     for (i = 2; i <= u[0]; i++) {
  1200.         if (u[i]) {
  1201.             v = int_new(u[0] -(i - 1));
  1202.             for (j = 1; j <= v[0]; j++) {
  1203.                 v[j] = u[i++];
  1204.             }
  1205.             int_free(u);
  1206.             return (v);
  1207.         }
  1208.     }
  1209.     /*  Return zero if all components zero */
  1210.     return int_con(0);
  1211. }
  1212.  
  1213. /* Not used */
  1214. int    value(char *s)                                                    /*;value*/
  1215. {
  1216.     /* Convert a numeric string to an integer. */
  1217.     return atoi(s);
  1218. }
  1219.  
  1220. static int    *int_new(int n)                                        /*;int_new*/
  1221. {
  1222.     /* Allocate new integer with n components, each initially zero */
  1223.     int       *p;
  1224.     int        i;
  1225.  
  1226.     p =(int *) ecalloct((unsigned) n + 1, sizeof(int), "int-new");
  1227.     p[0] = n;
  1228.     for (i = 1; i <= n; i++)
  1229.         p[i] = 0;
  1230.     return p;
  1231. }
  1232.  
  1233. static void int_free(int *n)                                    /*;int_free*/
  1234. {
  1235.     efreet((char *) n, "int-free");
  1236. }
  1237.  
  1238. int    *int_con(int i)                                                /*;int_con*/
  1239. {
  1240.     /* Build multi-precision integer from standard (short) integer */
  1241.     int       *p;
  1242.  
  1243.     if (i<-BAS || i > BAS) chaos("int_con arg out of range ");
  1244.  
  1245.     p = int_new(1);
  1246.     p[1] = i;
  1247.     return p;
  1248. }
  1249.  
  1250. int    *int_copy(int *u)                                        /*;int_copy*/
  1251. {
  1252.     /* Return copy of multi-precision integer u */
  1253.     int       *p;
  1254.     int        n, i;
  1255.  
  1256.     p = int_new(n = u[0]);
  1257.     for (i = 1; i <= n; i++)
  1258.         p[i] = u[i];
  1259.  
  1260.     return p;
  1261. }
  1262.  
  1263. int    int_eqz(int *u)                                                /*;int_eqz*/
  1264. {
  1265.     /*  Compare multi-precision integer with zero */
  1266.  
  1267.     int        i;
  1268.  
  1269.     for (i = 1; i <= u[0]; i++)
  1270.         if (u[i]) return (FALSE);
  1271.     return TRUE;
  1272. }
  1273.  
  1274. int    int_nez(int *u)                                                    /*;int_nez*/
  1275. {
  1276.     /* Compare multi-precision integer with zero; true if not zero */
  1277.     return ! int_eqz(u);
  1278. }
  1279.  
  1280. /* end module ada - arith */
  1281.  
  1282. #ifdef DEBUG
  1283. /* debugging and test procedures */
  1284. void int_print(int *u)                                            /*;int_print*/
  1285. {
  1286.     /* Dump multi-precision integer to standard output */
  1287.     int        i, n;
  1288.     n = u[0];
  1289.     if (n <= 0) {
  1290.         chaos("ill-formatted arbitrary precision integer - lengt <= 0");
  1291.     }
  1292.     printf("integer: %d components\n", u[0]);
  1293.     printf("%3d %*d\n", 1, DIGS+1, u[1]); /* allow for possible sign */
  1294.     for (i = 2; i <= u[0]; i++)
  1295.         printf("%3d  %0*d\n", i, DIGS, u[i]);
  1296. }
  1297. #endif
  1298.  
  1299. /*        module ada - arith   */
  1300. /*
  1301.  * All integers are stored in the multiple precision form used
  1302.  * in the package which is included as part of this program,
  1303.  * and rationals are stored as a pair [numerator, denominator],
  1304.  * with the denominator always positive, and [0] stored as [[0], n].
  1305.  * A rational plus or minus infinity may be represnted as
  1306.  * [n, [0]] or [-n, [0]] respectively. A rational of the form [[0], [0]]
  1307.  * will have an undefined effect if used in an operation.
  1308.  */
  1309.  
  1310. /* Macros for access to rational numbers: */
  1311.  
  1312. #define num(x) x->rnum
  1313. #define den(x) x->rden
  1314.  
  1315. Rational RAT_TWO;
  1316. Rational ADA_MAX_FLOAT;
  1317.  
  1318. /*
  1319.  * R A T I O N A L     A R I T H M E T I C     P A C K A G E
  1320.  
  1321.  *             Robert B. K. Dewar
  1322.  *              June 18th, 1980
  1323.  
  1324.  * This package contains a set of routines for performing
  1325.  * rational multiple precision arithmetic.
  1326.  
  1327.  *
  1328.  * A rational number is represented as a struct with two components, rnum
  1329.  * and rden, where num is the numerator, which may be negative, and den is
  1330.  * the denominator, which is non-zero and positive. Usually rational
  1331.  * numbers are reduced to lowest terms(see procedure rat_red), but none of
  1332.  * the routines depend on this assumption. The numerator and denominator
  1333.  * are represented as multiple precision integers, using the standard
  1334.  * formats for the multiple precision integer package, which must be
  1335.  * included as part of the program.
  1336.  * The macros num() and den() are used often, for comparison with the
  1337.  * original SETL code.
  1338.  *
  1339.  * The following routines are provided in the package
  1340.  
  1341.  *    rat_abs        rational absolute value
  1342.  *    rat_add        rational addition
  1343.  *    rat_div        rational division
  1344.  *    rat_eql        rational equal to
  1345.  *    rat_exp        raise rational to multiple precision integer power
  1346.  *    rat_fri        convert multiple precision integers to rational
  1347.  *    rat_frr        convert real to rational
  1348.  *    rat_frs        convert rational from string
  1349.  *    rat_geq        rational greater than or equal to
  1350.  *    rat_gtr        rational greater than
  1351.  *    rat_leq        rational less than or equal to
  1352.  *    rat_lss        rational less than
  1353.  *    rat_mul        rational multiplication
  1354.  *    rat_neq        rational not equal to
  1355.  *    rat_rec        rational reciprocal
  1356.  *    rat_red        reduce rational to lowest terms
  1357.  *    rat_str        convert integral fraction to string fraction
  1358.  *    rat_sub        rational subtraction
  1359.  *    rat_tor        convert rational to real
  1360.  *    rat_toi        convert rational to integer(rounds)
  1361.  *    rat_tos        convert rational to string
  1362.  *    rat_tup        convert string fraction to integral fraction
  1363.  *    rat_umin    rational unary minus
  1364.  */
  1365.  
  1366. /*
  1367.  * The following procedures are introduced in the translation to C:
  1368.  *    rat_new        allocate new rational, set components
  1369.  *    rat_free    free rational(deallocate space, including
  1370.  *              space used by components
  1371.  */
  1372.  
  1373. void rat_init()                                                    /*;rat_init*/
  1374. {
  1375.     /* Initialize rational and multi-precision packages */
  1376.  
  1377.     extern  Rational RAT_TWO;
  1378.     extern  Rational ADA_MAX_FLOAT;
  1379.  
  1380.     RAT_TWO = rat_new(int_con(2), int_con(1));
  1381.  
  1382.     /*    ADA_MAX_FLOAT    =
  1383.      *    [[79, 228, 162, 514, 264, 337, 593, 543, 950, 336], [1]],
  1384.      *                $ rational form of MAX_FLOAT
  1385.      *    ADA_MAX_INTEGER_MP= [1, 073, 741, 823],
  1386.      *                $ long form of ADA_MAX_INTEGER
  1387.      */
  1388.     /* Some C compilers do not like to see the most negative integer as
  1389.      * a constant, so we initialize ADA_MIN_INTEGER here by assingment
  1390.      * (assuming they can subtract properly!)
  1391.      */
  1392.     ADA_MAX_INTEGER = MAX_INTEGER; /* to be sure */
  1393.     ADA_MIN_INTEGER = -ADA_MAX_INTEGER;
  1394.     ADA_MIN_INTEGER--;
  1395.     ADA_MAX_INTEGER_MP = int_fri(ADA_MAX_INTEGER);
  1396.     ADA_MIN_INTEGER_MP = int_sub(int_umin(ADA_MAX_INTEGER_MP), int_fri(1));
  1397. #ifdef MAX_INTEGER_LONG
  1398.     MIN_INTEGER_LONG = - MAX_INTEGER_LONG;
  1399.     MIN_INTEGER_LONG--;
  1400.     MAX_INTEGER_LONG_MP = int_frs(MAX_INTEGER_LONG_STRING);
  1401.     MIN_INTEGER_LONG_MP = int_sub(int_umin(MAX_INTEGER_LONG_MP),
  1402.       int_fri(1));
  1403.     ADA_MIN_FIXED = MIN_INTEGER_LONG;
  1404.     ADA_MAX_FIXED = MAX_INTEGER_LONG;
  1405.     MAX_INTEGER_LONG_MP = int_frs(MAX_INTEGER_LONG_STRING);
  1406.     MIN_INTEGER_LONG_MP = int_sub(int_umin(MAX_INTEGER_LONG_MP),
  1407.       int_fri(1));
  1408.     ADA_MIN_FIXED_MP = MIN_INTEGER_LONG_MP;
  1409.     ADA_MAX_FIXED_MP = MAX_INTEGER_LONG_MP;
  1410. #endif
  1411. #ifndef MAX_INTEGER_LONG
  1412.     ADA_MIN_FIXED = ADA_MIN_INTEGER;
  1413.     ADA_MAX_FIXED = ADA_MAX_INTEGER;
  1414.     ADA_MIN_FIXED_MP = ADA_MIN_INTEGER_MP;
  1415.     ADA_MAX_FIXED_MP = ADA_MAX_INTEGER_MP;
  1416. #endif
  1417.     ADA_MAX_FLOAT = rat_new(
  1418.       int_frs("79228162514264337593543950336"),
  1419.       int_con(1));
  1420. }
  1421.  
  1422. Rational rat_new(int *u, int *v)                                /*;rat_new*/
  1423. {
  1424.     Rational r;
  1425.  
  1426.     r =(Rational) emalloct(sizeof(Rational_s), "rat-new");
  1427.     r -> rnum = u;
  1428.     r -> rden = v;
  1429.     return r;
  1430. }
  1431.  
  1432. static void rat_free(Rational r)                                /*;rat_free*/
  1433. {
  1434.     if (r->rnum != (int *)0) efreet((char *) r->rnum, "rat-free-num");
  1435.     if (r->rden != (int *)0) efreet((char *) r->rden, "rat-free-den");
  1436.     efreet((char *) r, "rat-free-rat");
  1437. }
  1438.  
  1439. #ifdef DEBUG
  1440. void rat_print(Rational r)                                    /*;rat_print*/
  1441. {
  1442.     printf("rational \n");
  1443.     int_print(r -> rnum);
  1444.     int_print(r -> rden);
  1445. }
  1446. #endif
  1447.  
  1448. Rational rat_abs(Rational u)                                        /*;rat_abs*/
  1449. {
  1450.     /* Absolute value of rational number */
  1451.  
  1452.     return rat_new(int_abs(num(u)), int_copy(den(u)));
  1453. }
  1454.  
  1455. Rational rat_add(Rational u, Rational v)                            /*;rat_add*/
  1456. {
  1457.     /* Add rational numbers
  1458.      *
  1459.      *   un        vn          un * vd  +  vn * ud
  1460.      *   --     +  --     =    -------------------
  1461.      *   ud        vd            ud * vd
  1462.      */
  1463.  
  1464.     int       *rn, *rm, *tn, *tm;
  1465.  
  1466.     tn = int_mul(num(u), den(v));
  1467.     tm = int_mul(num(v), den(u));
  1468.     rn = int_add(tn, tm);
  1469.     rm = int_mul(den(u), den(v));
  1470.     int_free(tn); 
  1471.     int_free(tm); /* free temporaries */
  1472.  
  1473.     return rat_new(rn, rm);
  1474. }
  1475.  
  1476.  
  1477. Rational rat_div(Rational u, Rational v)                            /*;rat_div*/
  1478. {
  1479.     /*
  1480.      * Divide rational numbers
  1481.      *
  1482.      *    un
  1483.      *    --
  1484.      *    ud
  1485.      *          un * vd
  1486.      *   ----    =      -------
  1487.      *          vn * ud
  1488.      *    vn
  1489.      *    --
  1490.      *    vd
  1491.      *
  1492.      * Test for division by zero
  1493.      */
  1494.  
  1495.     int       *rn, *rm;
  1496.  
  1497.     /* Return NULL instead of OM */
  1498.     if (int_eqz(num(v))) return ((Rational)0);
  1499.  
  1500.     /* Divisor is non-zero */
  1501.  
  1502.     else {
  1503.         rn = int_mul(num(u), den(v));
  1504.         rm = int_mul(num(v), den(u));
  1505.         /*   if rm < 0 then
  1506.          *    rn := -rn;
  1507.          *    rm := abs rm;
  1508.          *   end if;
  1509.          */
  1510.         if (rm[1] < 0) {
  1511.             rn[1] = -rn[1];
  1512.             rm[1] = -rm[1];
  1513.         }
  1514.  
  1515.         return rat_new(rn, rm);
  1516.     }
  1517. }
  1518.  
  1519. int rat_eql(Rational u, Rational v)                                    /*;rat_eql*/
  1520. {
  1521.     /* Test rational numbers for equality */
  1522.     int    *rm, *rn, res;
  1523.  
  1524.     rn = int_mul(num(u), den(v));
  1525.     rm = int_mul(num(v), den(u));
  1526.     res = int_eql(rn, rm);
  1527.     int_free(rn); 
  1528.     int_free(rm); /* free temporaries */
  1529.     return res;
  1530. }
  1531.  
  1532. Rational rat_exp(Rational u, int *ea)        /*;rat_exp*/
  1533. {
  1534.     /* Raise rational number to multiple precision integer power
  1535.      *
  1536.      *  If e >= 0:           If e < 0:
  1537.      *
  1538.      *    un    un ** e         un             1         ud **(-e)
  1539.      *    -- ** e = -------         -- ** e = --------------- = ----------
  1540.      *    ud    ud ** e         ud          (un/ud) **(-e)     un **(-e)
  1541.      */
  1542.  
  1543.     int       *e;
  1544.  
  1545.     if (int_eqz(num(u)))
  1546.         return int_eqz(ea) ? rat_new(int_con(1), int_con(1))
  1547.           : rat_new(int_con(0), int_con(1));
  1548.  
  1549.     e = int_copy(ea);        /* preserve value semantics */
  1550.  
  1551.     if (e[1] < 0) {
  1552.         u = rat_rec(u);
  1553.         e[1] = -e[1];
  1554.     }
  1555.  
  1556.     return rat_new(int_exp(num(u), e), int_exp(den(u), e));
  1557. }
  1558.  
  1559. Rational rat_fri(int *u, int *v)        /*;rat_fri*/
  1560. {
  1561.     /* Convert multiple precision integers to a rational number. */
  1562.  
  1563.     return rat_new(int_copy(u), int_copy(v));
  1564. }
  1565.  
  1566. Rational rat_frr(double u)                                            /*;rat_frr*/
  1567. {
  1568.     /* Convert a SETL real to a rational number.
  1569.      *
  1570.      * converts a floating number u to a pair of integers [p, y] such that
  1571.      * u = 2.0**p * float y
  1572.      * Here y satisfies 2**(N-1) <= abs y < 2 ** N
  1573.      * where N is the number of mantissa bits.
  1574.      * This essentially separates the fraction and the exponent.
  1575.      */
  1576.     /* The present C version is a straightforward translation of the SETL
  1577.      * code. Still unanswered is the question of whether real values will
  1578.      * be represented in C using doubles or floats; for now we assume reals
  1579.      * as floats. A more efficient version using library function 'frexp'
  1580.      * may be possible.
  1581.      */
  1582.  
  1583.     int        sgn;
  1584.     float   ub, lb;
  1585.     int        p, y;
  1586.     if (u == 0.0)
  1587.         return rat_new(int_con(0), int_con(1));
  1588.  
  1589.     if (u < 0) {
  1590.         u = -u;
  1591.         sgn = -1;
  1592.     }
  1593.     else {
  1594.         sgn = 1;
  1595.     }
  1596.     ub = pow2(ADA_MANTISSA_SIZE);
  1597.     lb = pow2(ADA_MANTISSA_SIZE - 1);
  1598.  
  1599.     p =(int)(log((double) u) / log((double) 2.0)) - ADA_MANTISSA_SIZE;
  1600.     /* estimate the exponent */
  1601.     u = u * pow2(-p);
  1602.     /* and adjust number */
  1603.     while(u >= ub) {
  1604.         u /= 2.0;        /* scale down */
  1605.         p += 1;
  1606.     }
  1607.     while(u < lb) {
  1608.         u *= 2.0;        /* scale up; */
  1609.         p -= 1;
  1610.     }
  1611.     y = sgn *(int) u;
  1612.     return rat_mul(rat_exp(RAT_TWO, int_fri(p)),
  1613.       rat_new(int_fri(y), int_con(1)));
  1614. }
  1615.  
  1616. /* Not used */
  1617. Rational rat_frs(char *s)                                            /*;rat_frs*/
  1618. {
  1619.     /* Convert a string representing a decimal fraction to the corresponding
  1620.      * rational number.  The string consists of a digit string, optionally
  1621.      * containing a decimal point and preceded by an optional sign.     If an
  1622.      * erroneous string is passed as an argument, then OM is returned.
  1623.      *
  1624.      * First step is to find number of digits after decimal point
  1625.      */
  1626.  
  1627.     int        dp, n, frac;
  1628.     int       *i;
  1629.     char   *t;
  1630.  
  1631.     n = strlen(s);
  1632.  
  1633.     t = strrchr(s, '.');
  1634.     frac =(t == (char *)0) ? 0 :(t - s + 1);
  1635.     /* find position of decimal point */
  1636.     dp = 0;
  1637.     if (frac != 0) {
  1638.         dp = n - frac;
  1639.         t = emalloct((unsigned) n, "rat-frs");
  1640.         *t = '\0';        /* mark end of(initially) null string */
  1641.         if (frac > 1)
  1642.             strncat(t, s, frac - 1);
  1643.         strncat(t, s + frac, (n - frac));
  1644.     }
  1645.     else
  1646.         t = s;
  1647.     /* Then number is converted as integer, and result is obtained
  1648.      * by supplying the appropriate power of ten as the denominator.
  1649.      */
  1650.  
  1651.     i = int_frs(t);
  1652.  
  1653. #ifdef XDEFER
  1654.     /* defer this while putting salloc in place DS (2-22-86) */
  1655.     if (frac)
  1656.         efreet(t, "rat-frs");        /* return t if allocated here */
  1657. #endif
  1658.  
  1659.     if (i == (int *)0)
  1660.         return (Rational)0;
  1661.     else
  1662.         return rat_red(i, int_ptn(dp));
  1663. }
  1664.  
  1665. int rat_geq(Rational u, Rational v)                                    /*;rat_geq*/
  1666. {
  1667.     /* Compare rational numbers: return true if u >= v, false otherwise. */
  1668.  
  1669.     return ! rat_lss(u, v);
  1670. }
  1671.  
  1672. int rat_gtr(Rational u, Rational v)                                    /*;rat_gtr*/
  1673. {
  1674.     /* Compare rational numbers: return true if u > v, false otherwise. */
  1675.  
  1676.     return (rat_eql(u, v) ? 0 : rat_lss(v, u));
  1677. }
  1678.  
  1679. int rat_leq(Rational u, Rational v)                                    /*;rat_leq*/
  1680. {
  1681.     /* Compare rational numbers: return true if u <= v, false otherwise. */
  1682.  
  1683.     return (rat_eql(u, v) ? 1 : rat_lss(u, v));
  1684. }
  1685.  
  1686. int rat_lss(Rational u, Rational v)            /*;rat_lss*/
  1687. {
  1688.     /* Compare rational numbers: return true if u < v, false otherwise.
  1689.      *
  1690.      *   un        vn
  1691.      *   -- <  --    iff   un * vd  <  vn * ud
  1692.      *   ud        vd
  1693.      */
  1694.     int    *rn, *rd, res;
  1695.     rn = int_mul(num(u), den(v));
  1696.     rd = int_mul(num(v), den(u));
  1697.     res = int_lss(rn, rd);
  1698.     int_free(rn); 
  1699.     int_free(rd);
  1700.     return res;
  1701. }
  1702.  
  1703. Rational rat_mul(Rational u, Rational v)                            /*;rat_mul*/
  1704. {
  1705.     /*
  1706.      * Multiply rational numbers
  1707.      *
  1708.      *   un        vn          un * vn
  1709.      *   --     *  --     =    -------
  1710.      *   ud        vd          ud * vd
  1711.      */
  1712.  
  1713.     int       *rn, *rm;
  1714.  
  1715.     rn = int_mul(num(u), num(v));
  1716.     rm = int_mul(den(u), den(v));
  1717.  
  1718.     return rat_red(rn, rm);
  1719. }
  1720.  
  1721. int rat_neq(Rational u, Rational v)                                    /*;rat_neq*/
  1722. {
  1723.     /* Test rational numbers for inequality */
  1724.  
  1725.     return ! rat_eql(u, v);
  1726. }
  1727.  
  1728. Rational rat_rec(Rational u)                                        /*;rat_rec*/
  1729. {
  1730.     /* Find reciprocal of rational number(number should not be zero). */
  1731.  
  1732.     int       *un, *dn;
  1733.  
  1734.     un = int_copy(den(u));
  1735.     dn = int_copy(num(u));
  1736.     if (dn[1] < 0) {
  1737.         un[1] = -un[1];
  1738.         dn[1] = -dn[1];
  1739.     }
  1740.     return rat_new(un, dn);
  1741. }
  1742.  
  1743. Rational rat_red(int *un, int *ud)                                    /*;rat_red*/
  1744. {
  1745.     /* Form rational reduced to lowest terms
  1746.      *
  1747.      * This procedure is given as arguments the numerator and denominator
  1748.      * of a rational value(as multiple precision integers). It returns
  1749.      * the rational formed by reducing these values to lowest terms.
  1750.      *
  1751.      * First a special case: zero is reduced to [[0], [1]] .
  1752.      */
  1753.  
  1754.     int       *i, *j, *t, *gcd, usign, dsign, rsign;
  1755.  
  1756.     if (int_eqz(un)) {
  1757.         return rat_new(int_con(0), int_con(1));
  1758.  
  1759.         /* Another special case: plus or minus infinity is reduced to [[1], [0]]
  1760.          * or [[-1], [0]] .
  1761.          */
  1762.     }
  1763.     else {
  1764.         if (int_eqz(ud)) {
  1765.             return rat_new(int_con((un[1] < 0 ? -1 :((un)[1] > 0 ? 1 : 0))),
  1766.               int_con(0));
  1767.         }
  1768.         /* Else we must compute GCD, using Euclid's algorithm. */
  1769.  
  1770.         else {
  1771.             usign = SIGN(un[1]); 
  1772.             dsign = SIGN(ud[1]);
  1773.             rsign = usign == dsign ? 1 : -1;
  1774.             i = int_copy(un);
  1775.             i[1] = i[1] >= 0 ? i[1] : -i[1];
  1776.             j = int_copy(ud);
  1777.             j[1] = j[1] >= 0 ? j[1] : -j[1];
  1778.             if (int_gtr(j, i)) {
  1779.                 t = i;
  1780.                 i = j;
  1781.                 j = t;
  1782.             }
  1783.             /* Steps of Euclid's algorithm, at each step, i is greater than
  1784.              * or equal to j, i is replaced by j and j is replaced by the
  1785.              * remainder of dividing i by j. The algorithm terminates when
  1786.              * j is zero, with i being the greatest common divisor.
  1787.              */
  1788.  
  1789.             while(int_nez(j)) {
  1790.                 t = j;
  1791.                 j = int_rem(i, j);
  1792.                 i = t;
  1793.                 /* [i, j] := [j, int_rem(i, j)]; */
  1794.             }
  1795.  
  1796.             gcd = i;
  1797.  
  1798.             /* Now reduce the original to lowest terms using the computed GCD */
  1799.  
  1800.             i = int_quo(un, gcd);
  1801.             j = int_quo(ud, gcd);
  1802.             /* adjust signs if needed */
  1803.             if (i[1] < 0) i[1] = -i[1];
  1804.             if (j[1] < 0) j[1] = -j[1];
  1805.             if (rsign < 0) i[1] = -i[1];
  1806.             return  rat_new(i, j);
  1807.  
  1808.         }
  1809.     }
  1810. }
  1811.  
  1812. #ifdef TBSN
  1813. char **rat_str(int **q)                                                /*;rat_str*/
  1814. {
  1815.     /* Convert tuple of multiple precision integers to tuple of
  1816.      * strings.
  1817.      * We interpret the argument q as a pointer to an array of pointers
  1818.      * to multi-precision integers.
  1819.      */
  1820.  
  1821.     char *emalloct();
  1822.     char *p;
  1823.  
  1824.     p = ecalloct(2, sizeof(char *), "rat-str");
  1825.     p[0] = int_tos(*q[0]);
  1826.     p[1] = int_tos(*q[1]);
  1827.     return &p;
  1828. }
  1829. #endif
  1830.  
  1831. Rational rat_sub(Rational u, Rational v)                            /*;rat_sub*/
  1832. {
  1833.     /* Subtract rational numbers
  1834.      *
  1835.      *   un        vn          un * vd  -  vn * ud
  1836.      *   --     -  --     =    -------------------
  1837.      *   ud        vd            ud * vd
  1838.      */
  1839.  
  1840.     int       *rn, *rm, *tn, *tm;
  1841.  
  1842.     tn = int_mul(num(u), den(v));
  1843.     tm = int_mul(num(v), den(u));
  1844.     rn = int_sub(tn, tm);
  1845.     int_free(tn); 
  1846.     int_free(tm);
  1847.     rm = int_mul(den(u), den(v));
  1848.     return rat_new(rn, rm);
  1849. }
  1850.  
  1851. double rat_tor (Rational u, int n)                                    /*;rat_tor*/
  1852. {
  1853.     /* TBSL: need to review and test this code        ds 11 sep 84 */
  1854.     double    real1;
  1855.     int    sgn, numpow, denpow;
  1856.     int    *nu, *du;
  1857.     int    *ub, *lb, p, *lbnd;
  1858.     int    *iquo;
  1859.     long    ntmp;
  1860.     double    log_bas, rpow, res;
  1861.     /* Convert a multiple precision real to a SETL real, if possible.
  1862.      * Make copy and work with positives
  1863.      */
  1864.     u = rat_red(num(u), den(u));
  1865.     arith_overflow = 0;
  1866.     nu = num(u);
  1867.     sgn   = SIGN(nu[1]);
  1868.     nu[1] = ABS(nu[1]);
  1869.  
  1870.     /* Check for 0 */
  1871.  
  1872.     if (sgn == 0 ) { 
  1873.         rat_free(u);
  1874.         return 0.0;
  1875.     }
  1876.  
  1877.     /* Check for overflow */
  1878.  
  1879.     if (rat_gtr(u, ADA_MAX_FLOAT)) {
  1880.         arith_overflow = 1;
  1881.         return 0.0;
  1882.     }
  1883.  
  1884.     /* To find an accurate floating representation of a rational number,
  1885.      * we normalize it so that
  1886.      *
  1887.      *        2 ** (N - 1) <= (num/den) < 2 ** N
  1888.      *
  1889.      * where N is the size of the mantissa.
  1890.      * We can then perform a long division, and know that the mantissa is
  1891.      * accurate to 2 ** (-N) i.e. to the last bit.
  1892.      */
  1893.     real1 = BAS;
  1894.     log_bas = log(real1) / log(2.0);
  1895.     ub = int_frl((long) pow2(ADA_MANTISSA_SIZE));
  1896.     lb = int_frl((long) pow2(ADA_MANTISSA_SIZE - 1));
  1897.  
  1898.     nu = num(u);
  1899.     du = den(u);
  1900.     numpow = (int)((double)((nu[0]-1)) * log_bas + log((double)nu[1])/log(2.0));
  1901.     denpow = (int)((double)((du[0]-1)) * log_bas + log((double)du[1])/log(2.0));
  1902.  
  1903.     p = numpow - denpow - ADA_MANTISSA_SIZE;
  1904.     if (p < 0 ) {                /* -p > 0 */
  1905.         nu = int_mul(nu, int_exp(int_fri(2), int_fri(-p)));
  1906.     }
  1907.     else {                    /*  p >= 0  */
  1908.         du = int_mul(du, int_exp(int_fri(2), int_fri(p)));
  1909.     }
  1910.  
  1911.     while (int_geq(nu, int_mul(ub, du))) {
  1912.         du = int_add(du, du);
  1913.         p++;
  1914.     }
  1915.     lbnd = int_mul(lb, du);
  1916.     while (int_lss(nu, lbnd)) {
  1917.         nu = int_add(nu, nu);
  1918.         p--;
  1919.     }
  1920.     iquo = int_quo(nu, du);
  1921.     ntmp = int_tol(iquo);
  1922.     real1 = sgn * ntmp;
  1923.     rpow = pow2(p);
  1924.     res = real1 * rpow;
  1925.     return res;
  1926. }
  1927.  
  1928. int    rat_toi(Rational u)                                                /*;rat_toi*/
  1929. {
  1930.     /* Convert the rational number u to a SETL integer.  The number
  1931.      * u is rounded by adding rational 1/2 or -1/2.     The int_quo
  1932.      * procedure is used to obtain the result of dividing the
  1933.      * numerator by the denominator.  The resuling extended integer is
  1934.      * then converted to a SETL integer using int_toi.
  1935.      */
  1936.  
  1937.     int        t, s;
  1938.     Rational r;
  1939.  
  1940.     t = num(u)[1];
  1941.     s = t == 0 ? 0 : t > 0 ? 1 : -1;/* get sign of u */
  1942.     /* s := sign num(u)(1);        $ get sign of u */
  1943.     /*    add or subtract 1/2(or 0) */
  1944.     r = rat_add(u, rat_new(int_con(s), int_con(2)));
  1945.     return int_toi(int_quo(num(r), den(r)));
  1946.     /* get quotient and convert it */
  1947. }
  1948.  
  1949. long rat_tol(Rational u)                                        /*;rat_tol*/
  1950. {
  1951.     /* Convert the rational number u to a (long)integer.  The number
  1952.      * u is rounded by adding rational 1/2 or -1/2.     The int_quo
  1953.      * procedure is used to obtain the result of dividing the
  1954.      * numerator by the denominator.  The resuling extended integer is
  1955.      * then converted to a SETL integer using int_toi.
  1956.      */
  1957.  
  1958.     int        t, s;
  1959.     Rational r;
  1960.  
  1961.     t = num(u)[1];
  1962.     s = t == 0 ? 0 : t > 0 ? 1 : -1;/* get sign of u */
  1963.     /* s := sign num(u)(1);        $ get sign of u */
  1964.     /*    add or subtract 1/2(or 0) */
  1965.     r = rat_add(u, rat_new(int_con(s), int_con(2)));
  1966.     return int_tol(int_quo(num(r), den(r)));
  1967.     /* get quotient and convert it */
  1968. }
  1969.  
  1970. #ifdef TBSN
  1971.     proc rat_tos (u, n);        /*;rat_tos*/
  1972.     
  1973.     /*
  1974.      * Convert a rational number u to a decimal
  1975.      * string with n places after the decimal point. The result
  1976.      * is correctly rounded and always has the specified number
  1977.      * of digits after the decimal place (even if they are zero)
  1978.      *
  1979.      * First acquire the sign (and then work with positive numbers)
  1980.      */
  1981.     
  1982.     int *q, *r;
  1983.     
  1984.     sn :
  1985.     = '';
  1986.     if num(u)(1) < 0 then
  1987.     num(u)(1) :
  1988.     = abs num(u)(1);
  1989.     sn          :
  1990.     = '-';
  1991.     end if;
  1992.     
  1993.     /*
  1994.      * Form result by multiplying by power of ten corresponding
  1995.      * to number of decimal places and then doing division.
  1996.      */
  1997.     
  1998.     un :
  1999.     = int_mul (num(u), int_ptn (n));
  2000.     ud :
  2001.     = den(u);
  2002.     int_div(un, ud, &q, &r);
  2003.     if int_gtr (int_mul (r, [2]), ud) then
  2004.     q:
  2005.     =int_add(q, [1]);
  2006.     end if;
  2007.     
  2008.     /*
  2009.      * Return result by converting this integer to a string and
  2010.      * then supplying the decimal point at the appropriate position.
  2011.      */
  2012.     
  2013.     q :
  2014.     = int_tos (q);
  2015.     if #q <= n then
  2016.     q :
  2017.     = (n + 1 - #q) * '0' + q;
  2018.     end if;
  2019.     
  2020.     return sn + q(1..#q-n) + '.' + q(#q-n+1..);
  2021.     
  2022.     end proc rat_tos;
  2023. #endif
  2024.  
  2025. #ifdef TBSN
  2026. Rational *rat_tup(char **q)                                            /*;rat_tup*/
  2027. {
  2028.     /* Convert tuple of strings to tuple of multiple precision integers.
  2029.      * We interpret tuple of strings as pointer to array of pointers to
  2030.      * strings.
  2031.      */
  2032.  
  2033.     char *ecalloct();
  2034.     Rational *p;
  2035.  
  2036.     p = ecalloct(2, sizeof(Rational ), "rat-tup");
  2037.  
  2038.     p[0] = int_frs(*q[0]);
  2039.     p[1] = int_frs(*q[1]);
  2040.     return &p;
  2041. }
  2042. #endif
  2043.  
  2044. Rational rat_umin(Rational u)                                    /*;rat_umin*/
  2045. {
  2046.     /* Rational unary minus */
  2047.  
  2048.     return rat_new(int_umin(num(u)), int_copy(den(u)));
  2049. }
  2050.  
  2051. static double pow2(int p)                                             /*;pow2*/
  2052. {
  2053.     double running, result;
  2054.  
  2055.     result = 1.0;
  2056.  
  2057.     if (p < 0) {
  2058.         running = 0.5;
  2059.         p = -p;
  2060.     }
  2061.     else {
  2062.         running = 2.0;
  2063.     }
  2064.  
  2065.     while(p != 0) {
  2066.  
  2067.         /* If p is odd then multiply result by running. */
  2068.  
  2069.         if (p % 2 == 1)
  2070.             result = result * running;
  2071.  
  2072.         /* Square running. */
  2073.  
  2074.         running = running * running;
  2075.  
  2076.         p /= 2;
  2077.     }
  2078.  
  2079.     return result;
  2080. }
  2081.